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New methods are presented that allow straightforward application of complex nonlinearities to finite 
element based rotor dynamic analyses. The key features are 1) the methods can be implemented with 
existing finite element or dynamic simulation programs, 2) formulation is general for simple application to a 
wide range of problems, and 3) implementation is simplified because nonlinear aspects are separated from 
the linear part of the model. The new techniques are illustrated with examples of inertial nonlinearity and 
torquewhirl which can be important in rubbing turbomachinery. The sample analyses provide new 
understanding of these nonlinearities which are discussed. 

INTRODUCTION 

The modeling and simulation of rotordynamic systems with rubbing is generally avoided in the design 
verification of turbomachinery, but is receiving more attention because rubbing is becoming unavoidable in 
closp tolerance designs. Some understanding of general rub behavior has been developed over the years 
using simplistic Jeffcot models and various techniques. Ehrich [1] determined stability bounds for full 
annular rub of a Jeffcot model with a circular but flexibly supported stator using a limit cycle assumption and 
dynamic equilibrium. A number of analyses such as those by Childs [2] and Muszynska [3] have treated the 
Jeffcot model with piecewise linear transient analysis or linearized analysis based on harmonic balance 
methods. The results are enlightening, but cannot be extended to complex systems without significant 
difficulty. 

Much work has been recently done to apply harmonic balance methods to more general systems. The 
work of Lau, Cheung, and Wu [4] is most notable because it treats aperiodic response by including multiple 
time scales in the assumed harmonic components. Unfortunately, all harmonic balance and related schemes 
suffer the debilitating requirement that the response frequency content or some fundamental response ratio 
be known a priori. Thus, one cannot analyze a given system for general response. 

The other remaining option is nonlinear transient simulation. Flexibility in modeling is achieved at the 
expense of computing time. However, with recent developments in computing power and well known 
methods of modal synthesis and model reduction, nonlinear transient simulation of real rotor systems is now 
a reasonable approach. Kascak [5] and Kascak and Tomko [6] have done some of the few analyses 
documented of transient rubbing with more complex rotors. Padovan and Choy [7] have more recently 
addressed the Jeffcot model with transient rub simulation, but with the added complexity of turbine blades 
modeled by static stiffnesses. These analyses have interesting approaches, but their application to general 
models is not discussed. 

General approaches to rotordynamic nonlinear problems are lacking. It is the intent of this work to show 
how nonlinear effects important in rubbing can be included in analysis of completely general rotor/housing 
systems with minimal effort, using generally available software augmented by relatively simple user-defined 
modules. By way of examples, it will be shown that the new techniques can be used to model nonlinearities 
of high complexity without major reformulation of a problem by exploiting the separability of linear and 
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nonlinear parts. The result is highly productive nonlinear analysis that is modular in nature, and thus easier 
to perform for the average rotor design analyst. Method formulation and examples will address the specific 
effects of inertial nonlinearity and torque, which invariably occur in rubbing rotordynamic systems. 
Incorporating general rub laws with flexible rotor blades and disks will be addressed in a later paper. 
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Dynamical matrix 
Load weighting matrix 
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Linear system damping matrix 
Load torque coefficient (N m sec/rad) 
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Forces on body in Cartesian system (N) 
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Linear system stiffness matrix 
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Rotor stiffness matrix with steady torque T 
Rotor stiffness matrix, total 
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Linear system mass matrix 
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THE GENERAL FINITE ELEMENT BASED APPROACH 


Today, every turbomachinery design organization has available to it a finite element analysis software tool. 
Rotordynamic analysis should exploit this capability to the fullest. Rotor, housing, disk, and blade models 
can be constructed using powerful solids modelers. The finite element models can be reduced by linear 
modal analysis or substructuring. Such a reduced matrix model, in the form of linear and symmetric mass, 
damping, and stiffness matrices, is an excellent base for nonlinear simulation. 

To continue with general rotor analysis, one must have complex eigenanalysis and nonlinear transient 
capability, which is generally available. For some, the finite element program may have this capability, while 
for others, one may have a dynamic simulation program (ACSL, CSMP, or others) with these options. It is a 
simple matter to format the matrix data between a finite element program output and a simulation program. 

It is assumed that the linear model includes flexibility and mass distribution for rotor, housing, disks, blades 
or other subcomponents, and possibly linear symmetric bearing or seal models. The entire system is 
described by the linear equations of motion. 

Mq + Cq + Kq = f (1) 

where the subcomponents may or may not be connected. Disconnection results in block diagonal form of 
the matrices. The displacements q consist of independent sets of physical or modal degrees of freedom 
(DOFs) for each unconnected subcomponent. Since physical motions are preferred for applying 
nonlinearities, any modal coordinates may require conversion to the physical domain using user-prepared 
pseudo-modal transformations 

x = 0q (2) 

x = 0q ( 3 ) 

where x is the full physical degree of freedom set, a different size than q. <t> contains some identity matrix 
block diagonals for subcomponents described by physical degrees of freedom, and other block diagonals 
that are made up of mode shape coefficients. Matrix partitioning is specifically avoided here to maintain 
generality of the ordering of q. The physical forces are subject to the pseudo-modal transformation as well, 

f = 0 T F (4) 


The modal transformation can be implemented by the use of linear constraint equations, which are 
available in most finite element programs, or by manual coding within the nonlinear evaluation routines that 
will be used later. Note that only the physical DOFs to be used for evaluating or applying nonlinear loads 
need be treated. 

Several simplifications have been made in the foregoing to make preparation of the linear core model 
straightforward. The assumption of linearity and matrix symmetry along with an assumption of modes 
computed by undamped analysis ensures real and equal left and right eigenvectors, and thus the same real 
matrix 0 in all auxiliary equations. This technique is popular and has been used by Noah [8] among others. 

Given the linear core model, asymmetric linear effects can now be added directly to M, C, and K matrices. 
Examples are linear models of Alford forces, gyroscopics, and bearing, seal, or impeller cross-coupling. 

When modes are used, some physical asymmetric effects may require the modal transformation if they act at 
physical DOFs that depend on modes. For example, for stiffness, 

K' A = < * >T K a 0 ( 5 ) 

where the asymmetric fully physical matrix is K A , and the corresponding physical/modal counterpart is K A - 
K a can then be added directly to the original symmetric stiffness matrix. 
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In the interest of generality, the linear model can be cast in first order form if required for compatibility with 
available software such as dynamic simulation packages. This form is 


y= Ay+Bf 


where 




( 6 ) 


(7) 


and the auxiliary equations if modal coordinates are used are equations (2), (3) and (4). M" 1 will always exist 
as long as the DOFs of the linear model are independent and the problem is well posed. 


Arbitrary nonlinear forces can be applied via the physical loads F and can be calculated using any state 
variable in x, x, q, q or any externally defined variable such as time, rotational speed, temperature, or 
pressure. The states are directly available from integration of equation (1) or (6) provided by the analysis 
code used and application of equations (2) and (3). Loads returning to the linear core simulation may need 
conversion through equation (4). 


As a prelude to the development of specific nonlinear force formulations, consider that all linear modeling 
has been done in an inertial Cartesian coordinate system. Possible DOFs are translations X, Y, Z and right- 
handed rotations 0 X , 8y, ® z - 


INERTIAL NONLINEARITY 
The Moment Equations 

The general equations of motion of a massive rigid body in space are extremely nonlinear in form. Even for 
the typical case of a rotor disk, where we may assume axial inertia l zz = I, transverse inertia l xx = lyy = l 0 , and 
all products of inertia are zero, we obtain the moment equations 
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which are written using Euler angles precession #, nutation 8, and spin #, and the axes x-y-z follow the body 
except spin as shown in Figure 1 . It is desired to cast these equations into the fixed Cartesian system of the 
linear model. First, the assumption that nutation 8 is small (sinfl ~ 8, cos 8 — 1 ) will be made. This is a very 
good assumption for real rotor dynamics problems where shaft excursions are small, and affords great 
advantages. Given this assumption, the appropriate moment transformation is 
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The small nutation assumption allows the nutation to be vectorially related to the Cartesian rotations 8 x 
and 8y, since they are also small, and small angles can be vectorially combined. Thus, 


8 x = 8cos<p , and (10) 

Inserting equation (8) into (9), collecting terms, and substituting equations (10) and their time derivatives, 
the moment equations in fixed Cartesian coordinates are 
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where s is the total rotational speed, with the small nutation assumption 


S = <p + 0 


(12) 


On the right hand side of equation (11a), only the first term is linear, and should normally be included in 
the linear model as lumped transverse inertia. For constant speed, the second term is the familiar gyroscopic 
linear skew-symmetric coupling used by most rotordynamicists, the third term is zero, and the fourth term is 
nonlinear, depending on nutation angle and rate. The fourth term is only zero for circular precession orbits 
(8 = 0) or planar motion (8 = 0), and in general cannot be assumed small. Although the portion 88 x is quite 
small, the portion 8 0 could be quite large because both 8 and <t> are proportional to speed. The approximate 
order of the whole term is order 8 3 s 2 , and should be compared with other cross-coupled stiffness loads in the 
system model to determine significance when general orbits are expected. 

The right hand side of equation (1 1 b) includes the linear angular acceleration term Is which may already be 
included as lumped polar inertia in the linear model. The corresponding DOF already in the linear model, 0 Z , 
can thus be interpreted as the integral of total local shaft speed. The remaining terms in equation (11b) are 
clearly nonlinear. They are zero for circular orbits and planar motion and are not necessarily small for other 
orbits, as both terms are of the order 0* s 2 . 


The Force Equations 


The force equation for a rigid body can be written directly in the Cartesian system without resort to the 
Euler system, thus 
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(13) 


Clearly, these equations are linear and the mass should be lumped in the linear model. 


Implementation 

The nonlinear inertial effects discussed above could be important in rub analyses because rubs have a 
complicated effect on the shaft orbit (nutation) and local shaft speed due to the induced impact and torque. 
The nonlinear terms are relatively straightforward to implement given the equations previously presented and 
a few more easily derived from them. Figure 2 shows the procedure to be implemented as a user-defined 
"nonlinear inertia finite element”. Note that all required Euler angles and rates are computed from the 
Cartesian rotations using the small nutation assumption, and singularities are avoided with logic checks. 

Most importantly, since the nonlinear moments are to be imposed on the model via the load vector F, the 
nonlinear terms are flipped in sign from that in equations (11) as if they were moved to the left hand side. 

Any number of these “finite elements” could be used to model the inertial behavior of multiple turbine disks 
and impellers in the same system model by multiple calls to a subroutine that implements Figure 2. 
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Example 


To illustrate the effect of inertial nonlinearity, a simple five DOF model was constructed using the 
techniques described above. The linear core model is shown in Figure 3, and the matrix form is that of 
equation (1) where the DOFs are 

qT = [ X Y e x 0 y e z ) (14) 

Note the different X and Y stiffness to induce nutation rate and the offset disk to produce nutation angle 
under imbalance loading. Note also that the realistic assumption of l 0 = 0.51 will negate the last term in 
equation (11b). Two cases were run: 

1 . Linear, constant speed s of 5000 rad/sec, with constant imbalance load of 20 N (linear gyroscopics 
included). 

2. Full nonlinear, initial speed of 5000 rad/sec, constant imbalance load of 20 N modified during 
simulation to follow the new total rotation angle 8 Z . 

Secondary effects on the imbalance load due to slight speed change and disk tilt are neglected here so that 
they do not obscure the inertial nonlinearities. Damping was included to settle the orbit more quickly, as 
initial position of the shaft was centered. The speed of 5000 rad/sec is near the third critical speed, which 
involves mostly nutation. 

Figure 4 shows the orbit results for Cartesian translations and rotations. Note that the orbit for rotations is 
by definition a polar plot of nutation given the small nutation assumption. The nonlinearity slightly increased 
the translation, and increased nutation by a factor of about three. This is despite fairly heavy rotational 
damping indicated by the relatively uniform elliptical shape compared to the translation orbits. The net result 
is a threefold increase in bearing load at this speed due to nonlinearity. A check of the induced speed 
variation showed miniscule effect. Most of the nonlinear effect is directly from nutation, and total neglect of 
the s terms would have given similar results. 

Admittedly, the model was tuned (high speed, ellipticity, nutation) for a severe nonlinear condition, for the 
sake of illustration. The results do indicate, however, that the nonlinearity of nutation can be important even 
for very small nutation. The nonlinearity of variable speed can be better evaluated in rub simulations, to be 
presented in a later paper. 


INCLUDING SHAFT TORQUE 
Formulation and Implementation 

Suppose that the disk in Figure 1 has torque applied to it about the z axis such that it follows the z direction 
for all time. This is a reasonable assumption for rubbing or for fluid-induced drive or load torques in 
turbomachinery. Whenever a nutation angle develops, components of the torque occur along the Cartesian X 
and Y. Using equation (9) and a drive torque vector T pointing in the +z direction, the resulting fixed 
Cartesian moments are, for small nutation, 
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(15) 


Thus, following torque acts as a cross-coupled stiffness effect on transverse rotation. Constant torque can, 
in fact, be included in the system model as a linear cross-coupled stiffness. Since equation (15) is written for 
load on the system, the cross-coupled terms should be flipped in sign before adding to the system stiffness 
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matrix. If torque is variable as in rubbing, equation (15) can be used directly to update the system load vector 

F. 


Example 

To show that the small nutation assumption has no significant effect on results, the fully nonlinear model of 
Vance [9] was exercised using the new method above for constant load torque (linear cross-coupling). The 
nature of the Vance model assures circular orbit and constant speed, so all inertial nonlinearities vanish and 
gyroscopics are included as linear cross-coupled damping. Complex eigenanalysis was used to construct 
root loci of the torquewhirl system as functions of speed. Figure 5 shows a speed dependent root locus for a 
load torque coefficient of C[_ = +5.625 N m sec/rad (+50 in lb sec/rad). See reference [9] for all other 
parameters and definition of terms. The frequency and load torque level (C|_ times the speed) at the point 
where the locus crosses the imaginary axis is the stability boundary. 

Upon comparison with the exact solution via Euler’s equations, the stability boundary frequency divided by 
the speed was found to closely match (within a few percent) the exact whirling speed ratio f for nutation 
angles up to about 10 degrees. This was expected, because f (see Appendix A) does not depend on nutation 
angle if it is small. The comparison of stability boundary torque computed by both methods also produced a 
very good match. This correlation forces disagreement with the comment in reference [9] that the inequality 

<t> + ¥= 0cos 6 + Ift (16) 

“is central to an understanding of how torquewhirl is produced”, because the linearized model assumes cos# 
is exactly unity. Rather, it is clear from inspection of equation (15) that torquewhirl is produced by a cross- 
coupled rotational stiffness projected into the transverse plane. It is clearly dependent on the mode shape of 
the system, which determines the degree of projection achieved. Yim, Noah, and Vance [10] have recently 
come to the same conclusion. Another conclusion of reference [10], summarized in Table 1, is supported by 
the current linearized analysis as well. 

However, real turbomachinery has both load and drive torques at different shaft locations, so the net effect 
cannot be determined without specific analysis. The method introduced here makes such analyses 
particularly easy for constant torques at multiple locations on general rotors, without resort to nonlinear 
analysis. 


Buckling Effects 

Yim, Noah, and Vance [10] and Cohen and Porat [11] addressed the effect of shaft torque on the shaft 
stiffness due to a buckling load effect. This effect is separate from the projection effect treated above, but 
does lead to similar results in that the previously linear shaft stiffness is modified by symmetric and skew- 
symmetric changes that have nonlinear dependence on torque level. To model this effect, the entire stiffness 
matrix related to the rotor requires change as a function of torque. This can be done using the methods 
presented here, by including the shaft stiffness change and resultant delta internal forces and moments in the 
system load vector F. 

Since the stiffness change operations essentially require individual shaft element influence coefficients, 
they can best be done with a finite element program with an enhanced shaft finite element. For torque that 
can be assumed constant, the torque-modified stiffness can be directly calculated at the outset and 
incorporated into the system stiffness matrix, which would then be constant and linear but asymmetric. See 
[10] and [11]. Then substructuring (with care) if desired and all the steps outlined previously could be taken 
to perform nonlinear analysis. 

The case of transient torque requires full reformulation of the rotor stiffness matrix at each time step, which 
some finite element programs can perform (ANSYS, for example). As an alternative to full reformulation, one 
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could use a linear fit of the shaft influence coefficients about a nominal operating torque T 0 . Thus the total 
rotor stiffness could be written 


K TOj = K(To> + (T - T 0 ) K 

where K(T 0 ) embodies the nominal torque effect and is included in the linear model, and 
A d*(T)| 

K - JT I 


(17) 


(18) 


The forces on the system due to the delta variation of T with time can then be included in the system load 
vector F by the term 


-(T-T 0 )Kq (19) 

Kean be easily evaluated off-line by computing K(T) at two values of T around T 0 , and then performing the 
operation 

K(T 0 + AT + ) - K(T 0 - AT') (20) 

K = 

AT + + AT- 


The method described above is ideal for rubbing problems where rub torques are small compared to steady 
drive torque. 

In many realistic cases shafts operate far below their buckling limit, and the torque buckling effect may not 
be important compared to the projection effect. Since it is now clear how to separate the two effects, the 
question deserves further investigation. 


CONCLUSIONS 

1. Using the power of existing finite element and/or simulation software, quite general nonlinear effects 
can be applied to rotor/housing models of arbitrary linear complexity. 

2. Casting nonlinear rigid body equations of motion into Cartesian form facilitates the understanding of 
the nonlinearities and simplifies their full implementation in general models. The assumption of small 
nutation affords great advantages. 

3. Inertial nonlinearity can be important in problems with high speed and non-circular orbits, such as 
found in rubbing phenomena. 

4. Casting the effects of projected following torque into Cartesian form using small nutation greatly 
simplifies torque modeling. Constant torque can even be handled as a linear cross-coupled stiffness. 

5. Torques have two separate effects: projection of moments and torque buckling. The relative 
importance of each depends greatly on rotor design, specifically on mode shape and torque buckling 
margin, respectively. Torque buckling is amenable to relatively simple modeling in cases where 
transient torques are small relative to a steady value. 

6. The combined effect of the items discussed in rubbing, blade loss, and other rotordynamic situations 
requires further study. This author will continue investigations, including modeling general bladed 
disks with rubbing. 
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APPENDIX A 


Stability of an Idealized Rotor with Following Torque 

The model of Vance [9], viscous case, was analyzed using Euler’s equations with the following results. See [9] for 
full definitions. 

Nonsynchronous Whirling Speed Ratio: 
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where w s = <t> + b = shaft speed 

Stability Boundary Torque 
Tb = -Cq l 2 w s f cos5* 

Note that neither quantity depends on 8 * if 8* is small (small nutation). 
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TABLE 1. Following torque stability influence 
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FIGURE 1. Coordinates and Euler angles for a rigid body 
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FIGURE 2. Algorithm for inertial nonlinearity finite element 
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FIGURE 3. Simplified rotor model for inertial nonlinearity example 




FIGURE 4. Response orbits of the simplified rotor model 
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FIGURE 5. Root locus for linearized Vance model 
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